An Asymptotic-Preserving all-speed scheme for the 
Euler and Navier-Stokes equations 



Floraine Cordier 1,2,3 , Pierre Degond 2,3 , Anela Kumbaro 1 

1 CEA-Saclay DEN, DM2S, SFME, LETR F-91191 Gif-sur-Yvette, France. 
floraine.cordier@cea.fr, anela.kumbaro@cea.fr 
2 Universite de Toulouse; UPS, INSA, UT1, UTM ; Institut de Mathematiques de Toulouse 

F-31062 Toulouse, France. 
3 CNRS; Institut de Mathematiques de Toulouse UMR 5219 ; F-31062 Toulouse, France. 

pierre.degond@math.univ-toulouse.fr 

^ Abstract 

We present an Asymptotic-Preserving 'all-speed' scheme for the simulation of 
compressible flows valid at all Mach-numbers ranging from very small to order unity. 
The scheme is based on a semi-implicit discretization which treats the acoustic 
part implicitly and the convective and diffusive parts explicitly. This discretiza- 
tion, which is the key to the Asymptotic-Preserving property, provides a consistent 
approximation of both the hyperbolic compressible regime and the elliptic incom- 
es — pressible regime. The divergence- free condition on the velocity in the incompress- 
ible regime is respected, and an the pressure is computed via an elliptic equation 
resulting from a suitable combination of the momentum and energy equations. The 
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implicit treatment of the acoustic part allows the time-step to be independent of 
the Mach number. The scheme is conservative and applies to steady or unsteady 
flows and to general equations of state. One and Two-dimensional numerical results 
J> provide a validation of the Asymptotic-Preserving 'all-speed' properties. 
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1 Introduction 

The numerical simulation of fluid flows at all Mach numbers is an active field of research. 
The occurrence of low Mach number regions in a globally compressible flow may be caused 
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by the boundary or initial conditions (e.g. in a fluid at rest subject to a supersonic jet ), 
by the geometry of the problem (e.g. in a nozzle with a large variation of the section), or 
by the underlying Physics (e.g. in the case of phase changes). This occurrence gives rise 
to specific numerical issues which are discussed below. 

When the Mach number tends to zero, compressible flow equations converge to incom- 
pressible equations: the compressible Euler equations in the inviscid case (respectively the 
compressible Navier-Stokes equations in the viscous case) converge to the incompressible 
Euler equations (respectively incompressible Navier-Stokes equations). This convergence 
has been studied mathematically by Klainerman and Majda [281 [29] (See also [9] IT8 | HT1 150] 
for reviews and references). However, in numerical simulations, it is very difficult to shift 
from compressible flow equations to incompressible ones in the regions where the Mach- 
number becomes very small. Therefore, it is necessary to design numerical methods for 
compressible flows that can handle both the compressible regime (i.e. local Mach-number 
of order unity) and the incompressible one (i.e. very small local Mach-number). This is 
the purpose of 'All-Speed schemes'. 

In this work, we derive an All-Speed scheme using the Asymptotic-Preserving method- 
ology. The Asymptotic-Preserving (AP) property is defined as follows. Consider a contin- 
uous physical model Ai e which involves a perturbation parameter e (here, e is the scaled 
Mach-number and M 1 represents the compressible Euler or Navier-Stokes model) which 
can range from e = C(1)toe<Cl values. Let A4° the limit of M. £ when e — > (here 
M.° is the incompressible Euler or Navier-Stokes model). Let now Ai A be a numerical 
scheme which provides a consistent discretization of Ai L with discrete time and space 
steps (At, Ax) = A. The scheme A4 £ A is said to be Asymptotic-Preserving (AP) if its 
stability condition is independent of e and if its limit Ai A as e — > provides a consistent 
discretization of the continuous limit model Ai°. The AP property is illustrated by the 
commutative diagram of fig. [T] 



£ -> 

A4 e * M° 



A — > 



A — > 



M 



A e^> 



Figure 1: Asymptotic-Preserving (AP) property: the upper horizontal arrow translates 
the assumption that the continuous model Ai £ tends to the limit model Ai° when e — > 0. 
The left vertical arrow expresses that Ai A is a consistent discretization of Ai e when the 
discretization parameter A — > 0. The lower horizontal arrow indicates that the scheme 
M. A has a limit Ai° A when £ — > for fixed A. Finally, the right vertical arrow expresses 
the AP-property: it says that the limit scheme Ai° A is a consistent discretization of the 
limit model A4° when A — > 0. 
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The present scheme is derived following the AP methodology and targets the situation 
of mixed flows where part of the flow has local Mach-number of order unity and is in a 
compressible regime and part of the flow has very small local Mach-number and is in 
the incompressible regime. More precisely, our scheme meets the following requirements. 
It is AP, i.e. it is consistent with both the compressible and incompressible regimes. 
The divergence-free condition on the velocity in the incompressible regime is explicitly 
satisfied up to the order of the approximation. The CFL condition is independent of the 
Mach-number. Therefore, the time-step is not constrained to be inversely proportional 
to the sound speed like. We remind that classical explicit schemes require such a time- 
step constraint which is very detrimental to the scheme efficiency in the small Mach- 
number regime. The scheme is conservative and preserves the correct shock speeds in the 
compressible regime. At last, the scheme applies to a general equation of state and to 
steady as well as unsteady flows. 

The present work is the continuation of earlier work on the construction of Asymptotic- 
Preserving schemes for fluid equations in the small Mach-number limit. In [15J, a first- 
order AP scheme is derived for the isentropic Euler equations. A second order version of 
this scheme based on the Kurganov-Tadmor central scheme methodology is proposed in 
|52j . Here, we extend the work of [15] to the full Euler and Navier-Stokes equations, i.e. 
including an energy equations instead of the isentropic assumption. This addition involves 
more than a simple technical adaptation. Indeed, the scheme has to be strongly modified 
in the choice of the terms that require an implicit treatment. Some of these terms have 
to be shifted from the mass to the energy conservation equation. With the use of a real 
gas equation of state, the resulting pressure equation becomes nonlinear and requires a 
specific treatment. We also provide a second-order extension of the method based on the 
classical MUSCL methodology which can apply to a larger software framework than the 
central scheme methodology. The numerical results will show that the passage to second 
order is qualitatively necessary to achieve a good accuracy. We also mention (21] which 
relates to [T5] but provides an alternate way of reaching the AP-property. 

Understanding why compressible flow solvers perform so poorly in the low Mach- 
number regime has triggered a vast literature since the seminal work of Chorin [6] . Volpe 
[57] observed that the numerical error increases when the Mach-number is decreased, at a 
constant mesh and that the convergence rate deteriorates noticeably. Guillard and Viozat 
[23] observe that an upwind space discretization leads to pressure fluctuations of the order 
of the Mach number e while in the continuous case the pressure fluctuations are of order 
£ 2 . This difference originates from the upwinding terms, and more precisely from the 
eigenvalues of the Jacobian matrix whose order of magnitude is the sound velocity. The 
argument has been developed further in [T6] . 

Compressible codes also require an increasingly large Computational time as the in- 
compressible regime gets closer. Indeed, the CFL stability condition for an explicit scheme 
reads At < |A ^ X| , where At is the time-step, Ax the space step, and A max is the fastest 
characteristic wave and can be written A max = u ± c, u being the fluid velocity and c 
the sound velocity. In scaled variables (see below for details on the scaling), the Mach 
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number £ appears explicitly in the stability condition as follows: 



Ax 



Ax 



Ax 



(1.1) 



At < — 



■=7 = £ 



max 



max |u ± - 

1 £ 



max I eu ± c 



where the tildes denote scaled quantities and the sound speed is now written c/e where 
c = 0(1). The time-step is therefore roughly proportional to the Mach number £ and is 
dramatically reduced when £ is small. 

The design of specific schemes for the small Mach-number regime has consequently 
triggered an abundant literature, following various tracks. A first track consists in apply- 
ing preconditioning methodologies. These methods have been initiated by the 'artificial 
compressibility' technique of Chorin [5] and consist in multiplying the time-derivatives by 
a suitable matrix. They aim at modifying the eigenvalues of the compressible system in 
order to reduce the disparity between the acoustic and fluid wave speeds [5] 13? ] 138 ] 154] , 155] . 
However, problems due to Computational instabilities related to the structure of the eigen- 
vectors [ID] and to the fact that the divergence-free constraint on the velocity is not always 
respected need to be dealt with. In most cases, these methods only apply to steady-state 
computations, since the time derivatives are modified. For non-stationary flows, dual 
time-stepping techniques can be introduced [1] to recover time-accuracy. Working with 
the original compressible equations, [TH] construct a semi-implicit Roe-type solver by 
decomposing the Jacobian matrix into the fast and slow eigenvalues, the former being 
treated implicitly. In |45] . the proposed scheme includes an implicit predictor convective 
step, followed by a semi-implicit corrector step. 

A second track consist in focusing on the pressure equation. To this aim, a natural 
idea is to adapt classical incompressible schemes to the compressible case. The pressure- 
correction method SIMPLE [2T]|4"T] solves an elliptic pressure correction equation obtained 
via the mass conservation equation and the equation of state. In [44] . the elliptic pressure 
correction equation is obtained by introducing the pressure equation (derived from the 
energy equation) in the momentum equation. These methods respect the divergence 
constraint on the velocity but the formulation is not always conservative. In the ICE 
(Implicit Continuous Eulerian) method introduced by Harlow and Amsden and followers 
[UEH], a splitting method is introduced between the explicit convective part and implicit 
acoustic part. However, the ICE method is not conservative and inaccurate shock speeds 
are observed. Klein [30J proposes a semi-implicit scheme which solves explicitly the leading 
order contribution of the pressure and the lower orders, implicitly. Other ways generating 
elliptic equations on the pressure can be found in [331 1131 EE SHI ESI EH] ■ 

A third track consists in using gauge (or Hodge) decomposition of the flow variables [7]. 
Indeed, the incompressible velocity between divergence free, it is tempting to decompose 
the compressible velocity into a divergence-free and a curl-free part. Semi-implicit time 
discretizations are used for the divergence-free part. The gauge decomposition was used in 
an earlier attempt to derive and AP-scheme [13]. However, the method was too complex 
and never used. 

To some extent, our work belongs to the second class and relies on the introduction 
of a suitable elliptic equation on the pressure. However, it departs from previous work 
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in that the problem is discretized in a single step, which reduces the Computational cost 
compared to predictor-corrector procedures, that the scheme is conservative and that the 
only equation solved implicitly is the elliptic equation, whose construction is extremely 
simple. 

More generally, AP-schemes have previously been proposed for neutron transport prob- 
lems [34J, multiscale kinetic equations [21], hyperbolic heat equations [22J, relaxation limit 
of hyperbolic models [40], plasmas in the quasi-neutral limit [HIE] or in the large magnetic 
field limit [12]. 

The outline of this paper is as follows. We first provide a semi-implicit AP time 
discretization of the compressible flow equations in section [2} Then, we derive the fully 
discrete (in time and space) AP-scheme at first order in section [3] The construction of 
an elliptic equation on the pressure as well as the resolution of the scheme is detailed, 
and the outline of the extension to a second order scheme is given. Then, we perform 
the asymptotic analysis of the proposed scheme in section |4j in order to show the AP 
property. Numerical results presented in section [5] provide a validation of the scheme in 
both the compressible and close-to-incompressible regimes. Finally, a conclusion is drawn 
at section [6j 



2 Time semi-discrete scheme 



We start with the Navier-Stokes equations (2.1)-(2.3): 
3 t p + V • pu = 0, 

a t pu + v ■ (pu <g> u) + vp = v • 



3 tP E + V- (pHu) 



1 



pv 
A 

Cr 



(Vu + Vu T )--(V-u)I 



+ pf< 



cxt ) 



-VH 



+ pf« 



cxt 



W= pE = -pu z + ph-p, 



(2.1) 
(2.2) 

(2.3) 

(2.4) 



where p is the density, u is the velocity, p is the pressure, h. is the enthalpy, E is the total 
energy, H = E + Ms the total enthalpy, v is the kinematic viscosity, Vu T is the transpose 
of the gradient of the velocity, I is the identity matrix, A is the conductivity, C p is the 
specific heat capacity, and f ext represent external forces like gravity. The contribution of 
the term V- [pv ((Vu + Vu T ) — |(V • u)l)] u in the energy equation has been neglected, 
according to the models used in the CEA codes FLICA4 0, [53] and CATHARE [42], but 
could easily be added. We consider a general equation of state linking the density, the 
pressure and the enthalpy: 

p = p(p,H). (2.5) 

In this paper we deal with scaled equations. The scaling parameters po, po, u-o, *o are 
introduced along with the scaled variables, denoted by a tilde. 



P_ 
Po ; 



u 



u 

u " 



V_ 
Vo' 



X 

Xo' 



Po 
Po 



E, H 



Po 1 
Po 



(2.6) 
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The scaled equation are the following (we will omit the tildes in the remainder of the 
paper): 



3 t p + V ■ pu 
9 t pu + V ■ (pu <g> u) + -^Vp 

£ z 

3 tP E + V- (pHu) 
W= pE 



o, 



Re 



1 



Re - Pr 

-e 2 pu 2 + ph. 



Vu + Vu T ) - 
A h + £ 2 pf ext • u, 
P> 



2 



(V ■ u)I 



+ pf< 



cxt y 



where the parameters resulting from the scaling are: 

u x 



2 _ Po u o 
Po 



Re = 



Pr = 



Po^C. 



(2.7) 
(2.8) 

(2.9) 
(2.10) 

(2.11) 



The parameter £ represents a global Mach number characterizing the flow and the nondi- 
mensionalisation. It is different from the local Mach number. The parameter Re is the 
Reynolds number and Pr is the Prandtl number. 

For the sake of simplicity, the scheme is presented on the full Euler equations, which 
represent the convective part of the Navier-Stokes equations. The right-hand terms in the 
Navier-Stokes equations (2.1 )-(2.3) will be included later in explicit source terms, and the 
time semi-discretization will not be modified. 

The AP time semi-discrete scheme is written as follows: 



^n+l 



q 



At 

n+1 



JL- + V ■ q n = 0, 



At 

W n+1 - 



V 



P n 



W 1 



At 



+ V ■ H n q n+1 
1 



+ ocp n ) + 
= 0, 



1 



OC£ 



VP 



n+l 



p n+1 e n+1 + -£ 2 p n (u n ) 2 



,n+1 H n+1 



p n +1 + ^ £ 2 p u (u n )2) 



(2.12) 
(2.13) 

(2.14) 
(2.15) 



where At is the time-step, t n = nAt and the superscript 'n' denotes the approximation 
of the variables at t n , q = pu, W = pE, and a is a small ad-hoc parameter independent 
of the Mach number and such that a ~ 1 , designed to prevent spurious oscillations in 
strong shock cases [IS]. The time discretization of the total energy W = pE splits into 
an implicit evaluation of the internal energy pe = ph — p, and in an explicit evaluation 
of the kinetic energy. The discretization of the space derivatives is detailed in the next 
section (Section pi). 



Let us make a few comments on the proposed scheme. First, the scheme being con- 
servative, we expect good shock properties in the compressible regime. Then, we will see 
that the implicit treatment of the pressure in the momentum equation (2.13) is a key 
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to the asymptotic preserving property (Section [42). An other noticeable feature is the 
implicit treatment of the momentum q in the energy equation, allowing us to construct 
an elliptic equation on the pressure. We now detail the resolution of the scheme and the 
construction of this elliptic equation. 

The scheme can be solved through the following steps: 

First, the density p n+1 is obtained via the resolution of the explicit continuity equation 



(2.12) 



An elliptic equation on the pressure is then solved. To construct this equation, the 



momentum equation (2.13) is rewritten as: 



q n+1 = q n - At V 



. q n <g) q n 



P n 



+ ocp n ) - At 



1 — at 1 



VP 



n+l 



This expression is inserted into the energy equation (2.14) and leads to: 

1 — <X£ 2 



W n+1 - At 2 



y(H n VP n+1 ) =<Mp 1 \q n ,W 1 ), 



(2.16) 



(2.17) 



where the right hand side § is explicit and is equal to: 



4)(p n , q n , W 1 ) =W n - AtV ■ H n q n + At 2 V ■ H n y • 



. q n <g> q n 



+ ocp r 



(2.18) 



Two cases can be considered : the specific case of a perfect gas equation of state, and 
the case of a general equation of state (EOS). 



Perfect gas EOS For a p erfect gas of polytropic constant y, the internal energy reads 
pe = ^jj-p. We can rewrite (2.17) as follows: 

1 - cxe 2 



p n+i _ [y _ 1)At 2 _ v . ( H - y p- +1 ) = *(p-,q n , W*), 



(2.19) 



with 



*(p n , q*, W*) = (y - 1 )4)(p n , q n , VT) - \{y - 1 )e 2 p n (u n ) 2 . (2.20) 



Equation (2.19) is an elliptic equation on the pressure. It allows us to find the pressure 



) n+1 , and then W n+1 . 



General EOS For a general equation of state, the internal energy reads pe = ph. — p. 
In this case, the following system has to be solved: 



,n+1 h n+1 



1 — CCF ^ ~ 

p n+i _ At 2 _ y . (H n y p n+1 ) = $'{p n , q n , W 11 ) 



(2.21) 



n+l v.n+1- 



P(p n+ ',H 



,n+l 
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where 



4>'(p n ,q n , 



-e 2 p n (u n ) 2 



AtV ■ H n q n + At 2 V • H n v 



q IL <g> q ,: 



+ <xp T 



(2.22) 

This still leads to an elliptic equation for the pressure, and the enthalpy is constrained by 
the value p n+1 of the density found by the resolution of the explicit continuity equation. 
Solving this system allows us to find p n+1 , H n+1 and W n+1 . 

The momentum q n+1 is finall y obt ained via the momentum equation (2.16), as p n+1 



is 



now known. Let us note that in ( 2.16[) a ll terms are 0(1 ). Indeed, we have 1 ^ VP n+1 



0(1) due to the elliptic equation (2.21) which implies that p n+1 = 0(e 2 ) in the Sobolev 
space H 2 given the elliptic regularity theorem, and using appropriate boundary conditions. 
We thus get that Vp n+1 = 0(e 2 ) . 

The proposed scheme presents two notable differences with the scheme for the isen- 
tropic equations presented in [15]. First, the density is taken explicitly in the continuity 
equation. Then, the elliptic equation is obtained by the insertion of the momentum equa- 
tion into the energy equation instead of into the continuity equation in the isentropic case. 
This difference is a consequence of the asymptotic analysis of the continuous full Euler 



equations (Section 4.1) where the divergence constraint on the velocity in the low Mach 
number regime is obtained from the energy equation. 



3 Full time and space discretization 

We present the full time and space discretization of the scheme for a first order scheme 
in a first part. Then we will extend the discretization to a second order scheme. We also 
insist on the centered space discretization of the implicit pressure. 
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3.1 First order scheme 



In the finite volume framework, the first order space discretization of the scheme for a 
general multidimensional system on a structured or unstructured mesh is given by: 



pr 1 - p? 

At 



qr +1 - q? 
At 



HI; 



III; 



c£ + q n 



^ + D n 

2 iv p 



S n 



veu(ij 

p n+1 h n+l _ p n+1 _ ^ 

At 



2e 2 



(3.1) 
(3.2) 



m, 



veu(i) 



i_in„n i T_in„n H n S n -I- H n S n 



III! 



At^H 



^ + Vl w + At- 

1 — CX£ 2 



m, 



uev(v) 



tv+1 i „n+1 • 



2e 2 



(p v n+ '+P 



(3.3) 



where, to simplify, we have introduced the notations: 



3[ 



c v n i; 



Vi 1 

1 qr®q? qy®q? . 

2 l pr p? 



+ <* z ^Div a > 



(3.4) 
(3.5) 



v[i) is the set of neighbors of the cell i, ni v is the unitary normal of the face between 
the i and v cells, is the surface of this face, Vj. is the volume of the cell i, and c{ v is 
+ 1 for an incoming normal of the face iv into the cell i and —1 for an outgoing normal, 
D£, = (D£, p , D£, , D£, w ) is the upwinding between the i and v cells, taken at the time n, 
and detailed below. 



system (2.12)-(2.14). 



Note that the energy equation (2.14) has been replaced by a discretization of the 



elliptic equation (2.21) on the pressure, the system so constituted being equivalent to the 



General source terms S n have been added and can include external forces such as 
gravity and the diffusive terms of the Navier-Stokes equations. 



Upwinding Centering the pressure term 1- y 2 Vp n+1 in the spatial discretization is a 
crucial feature of the low Mach number scheme. It does not affect the stability as it 
is an implicit term. Then, the upwinding only concerns the explicit part of the flux in 



the equations (3.1)-(3.3) and the eigenvalues of the Jacobian matrix of the corresponding 
system are: 

aa^ , |uj , u n + a/ aa 2 ra , (3.6) 



Un 
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where = u • n and a m is the sound speed defined by: 



3p p 3h 

The CFL condition for the stability of the scheme is: 



«-=ll/j£ + r{£l ■ (3-7) 



Ax 

At < (3.8) 

max(u n ± a/ aa 2 m ) 

Therefore the time-step At does not depend on the Mach number e contrary to a standard 
explicit method, as explained in the introduction. The time-step is based on the fluid 
velocity only : it does not take into account the acoustic velocity, which tends to infinity 
when the Mach number tends to zero and is responsible for the dramatic decrease of the 
time-step in the low Mach number regime. Also, the inaccuracy of explicit upwinding 
schemes is caused by the upwinding being based on the acoustic velocity, as recalled in 
the introduction and detailed in [23]. To avoid introducing wrong pressure fluctuations, 
we must keep the parameter cc small compared to \. 



In our scheme, a Lax-Friedrich upwinding is used. The term D| V in the discretization 



(3.1)-(3.3) gives the upwinding between the cells i and v and its expression is: 

D iv • n iv = I D iVq I • n lv = --(A™ ax UV v - V t ), (3.9) 

where V = (p,q, W) is the vector of conservative variables and 

(A™ x ) iv = max ( |u a | i +^4 i , |uj v + )• (3.10) 

Resolution of the discrete system Let us detail the steps in the resolution of the 
scheme. 



First, the mass equation (3.1) can be solved explicitly, and p n+1 is obtained. 



Then we solve the elliptic equation (3.3). For a perfect gas EOS, this elliptic equation 



is a linear system on the pressure and can be solved by inverting the system. In the case 



of a general EOS, the system constituted by the elliptic equation (3.3) and the equation 
of state is solved by means of a Newton method where the unknowns are the pressure and 
the enthalpy. 

We will note (p (q) , h (q) ) the pressure and enthalpy found by the q th iteration in the 
Newton method in order to find (p n+1 , h n+1 ) at the time t n+1 . Two iterations q and q + 1 
of the Newton method are linked by the following relation: 

= ($!) -fV q w q) rMp (q un on) 
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where the algorithm is initialized with (p'°',fi' ') = (p n ,h. n ). The function f is a vector 
defined as: 

f(p ( <>, h^) = [E 2 W q \ri q) )Mv [,, \rt q) )), (3-12) 



The first component fi comes from the elliptic equation (3.3) and the second component 
f 2 expresses the condition over the pressure and the enthalpy given by p n+1 : 



f 2 (p (q \h (q) 



,u+1 



P(p tq) ,h [ 



(3.13) 



In practice, the first component of f is £ 2 fi in order to avoid the division by the small 
parameter e in the term 2 " 2 £ . 



The matrix f (p q ,H q ) in (3.11) is the following: 



,23^ 2 3fi 

3p 3p 
3p 3h 

,n+1 „„ j AA/n+1 



(3.14) 



Solving the elliptic equation allows us to find p , h. and W n+ . Finally, the 
momentum equation (3.2) is solved to obtain q n+1 and u n . 



3.2 Second order scheme 

The first order scheme being too diffusive, we propose a second-order space discretization 
of the scheme. 



In the first order system, the full time and space discretization (3.1)-(3.3) could be 
written as: 

V n+1 — V n 

— l - + Y_ 0(Vf' n+1 , V^' n+1 ) = 0, (3.15) 

veu(i) 

where O is the numerical flux and Vi is the vector of conservative variables in the center 
of the cell i. The second order space discretization consists in evaluating the numerical 
flux (I> in the reconstructed and limited states V[ v and V? , which correspond to the 
vectors of conservative variables Vi and V v on the face between the cells i and v. We 
thus replace ®(VT l ,V™) by O ^(V[ v ) n , (V? ) n j . The minmod limiter is used to avoid 
spurious oscillations. 

On a two-dimensional Cartesian mesh, the reconstructed and limited states V L , . and 
V R . are given by the following expressions ([39], I2T]): 



y\ + y = Vij + ^minmod(V g - V w j, V i+ljj - Vy ) , (3.16) 
Vf +1 i = V l+1)j - lminmod(V 1+2)j - V l+1J , V i+1J - V t , j) , (3.17) 



where the minmod function is: 



minmod(x,-y) = -[sign(x) + sign(y)]min(|x|, |y|). (3.18) 
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The upwinding of the scheme is now given by: 



D 
D 
D 



D iv • n iv = ( D iVq | • n iv = -^(A^iv (V? v - V\ v ) , (3.19) 



IV,, 



where (A™ ax )i V is still given by equation (3.10). 



To solve this scheme, the reconstructed and limited states (V^,) n and (V^) n are first 
calculated from the conservative vector V n . In addition, we also obtain the corresponding 
pressure. These values are used to evaluate the numerical fluxes, the upwinding and the 
source terms. The mass equation is first solved explicitly and p n+1 is found. Then, the 
elliptic equation is solved by means of a Newton method, as explained for the first-order 
scheme. The momentum equation is then solved and V n+1 is obtained. 

A second-order discretization in time was intended via a Runge-Kutta method com- 
bined with Crank-Nicolson (RK2CN). However the scheme tended to generate spurious 
oscillations if the pressure was fully implicit (cx = 0). To date, the problem has not been 
solved by the authors. As the full implicit treatment of the pressure is necessary in the 
low Mach regime, the time and space second-order scheme could not be used. In the next 
part of the article, the "second-order scheme" will thus refer to the first-order in time and 
second-order in space scheme. 



4 Asymptotic preserving property 

Let us now show that the proposed scheme is asymptotic preserving. The asymptotic 
preserving property has been defined in the introduction. We first recall the asymptotic 
study of the full Euler equations as the methodology is used in the study of the asymptotic 
preserving property of the scheme. 



4.1 Asymptotic analysis of the continuous Euler equations 

Let us now investigate the limit of the full Euler equations when e — > 0. The method 
differs from the isentropic case [15] as the condition on the divergence of the velocity in 
the low Mach number regime is obtained via the energy equation instead of the continuity 
equation. This is a consequence of the density depending both on pressure and enthalpy 



(2.5). The analysis below extends the asymptotic analysis led by Klein in (301 |31] for the 



full Euler equations to a general equation of state. 

If we write the expansions of the variables p, p, u, H and W in powers of the Mach 

number e, e.g. p = p + £P(i) + £ 2 P(2) + • . ., and insert them in the full Euler equations, 
the leading order equations are: 

3 t P(0) + V • (p(0)U ( o)) = 0, (4.1) 

Vp (0 ) = 0, (4.2) 

3tW (0) + V • (p ( o)H (0) u ( o)) = 0, (4.3) 



12 



and the second order equation for the momentum is: 

9 t (pu) (0 ) + V(p(o)U (0 ) <8> u ( o)) + Vp (2 ) = 0. 



(4.4) 



The variable pp] is a dynamic pressure as it is directly linked to the speed of the fluid, 
while p(o) is a thermodynamic pressure as it appears in the energy equation. Eq.(4.2) 



yields that p(o} is independent of space. We assume that the boundary conditions are 
chosen such that the constant p(o) is independent of time. As the parameter e 2 appears 



in the expression (2.10) of W, at leading order we have: 

W (0 ) = P(o)e(0) and P(o)H (0 ) = P(o)e(o) +P(0), ie H (0 ) = H (0 ). 



(4.5) 



We first recall the simpler case of a perfect gas, then extend the analysis to a general 
equation of state. 



Perfect gas case For a perfect gas with a constant y, we have pe 



p. There- 



fore W, 



(0) 



rrPco] and P(0)H( ) 



Y-1 



P(0) are independent of space due to (4.2), and 



independent of time. The leading order of the energy equation (4.3) gives the divergence 
condition on the velocity in the zero Mach number limit: 



V • U( 0) = 



(4.6) 



General EOS case We drop the subscript (0) for simplicity. With W = ph— p, H = h 

(4.7) 



and (4.2), we get: 



3 t W+ V ■ (phu) = (3 t + u • V)(ph) + ph(V • u) = 0. 
Now, ph = ph(p,p) for a general EOS, and we get 

(9 t + u • V)(ph) = ^(3 t + u • V)p + ^(3 t + u • V)p 
3p 3p 

= -p^-V-u, 
3p 



(4.8) 
(4.9) 



thanks to ( |4.2[ ) and the assumption that 3 t p = 0. We collect the above equations and get 

(4.10) 



(ph-p^)V-u = 0, 
3p 



3ph 



,2 3h 



or, since ph - p^ = -p' 3p , 



3h 
3p 



Vu = 0. 



With the assumption that |^ ^ 0, we get the incompressibility condition: 



V-u = 0. 



(4.11; 



(4.12) 
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The divergence of the velocity being zero, the mass equation (4.1 ) becomes 



9 t p + u- yp = 0, 



(4.13) 



which expresses that the density is constant along a trajectory of any fluid element. By 
contrast, in the isentropic case, the low Mach number limit leads to a constant density in 
space. 



Equations (4.2), (4.4), (4.12) and (4.13) form the incompressible limit of the Euler 

then Metivier and Schochet in EH 



equations. Klainerman and Majda in [28, 29] , then Metivier and Schochet in [41] have 
shown that the solution of the compressible Euler equations converges towards the solution 
of the incompressible Euler equations when e tends to zero. 



4.2 Study of the asymptotic preserving property of the scheme 



Let us now show that the proposed scheme is asymptotic preserving. We expose the 
reasoning on the time semi-discrete scheme (2.12)-(2.14) for the sake of simplicity and 
readability. The extension to the full time and space discretization is straightforward. 
To show the asymptotic preserving property, we have to write the limit discrete scheme 



M° A when e 



and show that it is consistent with the continuous limit model M. at 



£ = 0. 

The continuous limit model Ai° is the following: 

9tP(0) + V ■ (P(o)U(o)) = 0, 
Vp(o) = 0, 

3 t (pu)(0) + V ■ (P(0]U(0) ® U( j) + V7T = 

9 t W( ) + V ■ (p ( o)H(o)U ( o)) = 0, 
H(o) = H(o), W(o) = P(o)6(o) = P(o}b(o) 



(4.14) 



P(0)J 



where n is a dynamic pressure and p(oj a thermodynamic pressure, and under the assump- 
tion that p(o) is independent of time. 

We introduce the expansions in powers of £ in the semi-discrete scheme (2.12)-(2.14) 
in the same way as in the asymptotic analysis of the continuous case (Section 4.1). Con- 
sidering the leading order equations and the equation of order two for the momentum 
equation, we obtain the discrete limit system A4° A : 



,n+1 
J (0) 



'(C) 



At 

(Pu)fo! 1 



+ V-(pf 0) uf 0) )=0, 



(01 



At 

(pe)-; 1 - (peft 



+ V- (pfoufo <g> u 



(0) 



(0) 



At 



n+1 u n+1 ' 



>P(0)6(0]. 



+ V ■ Ot-(o)P(o] "(o) 

n+1 f_ -in+1 



o, 



(4.15) 



n+1 
P(0) • 
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System (4.15) is clearly consistent with system (4.14). Therefore, the scheme is asymp 



totic preserving. 

Nonetheless, we show directly that (4.15) is also consistent with the incompressibility 
constraint, namely that 



Proposition: V ■ uT^t 1 = O(At) , where O(At) is independent of £. 



Remark 1: From now on, we drop the subscript (0) and O(At) will denote terms 
estimated by CAt with C independent of e. 



Remark 2: From (4.15), we deduce in particular that 



u 



,n+1 = p n +0 (At), 



n+1 



u n + 0(At), 



( P h) 



n+1 



(ph) n + p n+1 -p n + 0(At), 



with O(At) independent of £. From (4.17), we deduce that 

V-u n+1 = V-u n + 0(At). 



(4.16) 
(4.17) 
(4.18) 



(4.19) 



However, even if V - u = 0, this does not prove that V ■ u n+1 = O(At), since summing 
over all time steps will lead to V • u n+1 = 0(1 ). Therefore, we need to show directly that 
V ■ u n+1 = O(At). The proof is similar as in the continuous case. 



Remark 3: From the second equation of (4.15), we deduce that p n+ is independent of 
x. We assume that the boundary conditions are such that p n+1 is also independent of n, 



i.e. p n+1 = p n 



Proof: 



We write the fourth equation of (4.15) as 

( P h.r +1 -(phy 



At 



p n+l _ p n 
At 



+ u 



n+1 



V(H n p n+l ) + h n p n+1 V-u 



.n+1 



0. 



Since 



,n+l 



p n = 0, ( ph) n+1 = ( p h) { p n+1 , p n+1 ; 



(phr = (ph)( P n ,p n ), 



(4.20) 



(4.21) 
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we have, using (4.17) and the first equation of (4.15): 



(ph) n+1 - (ph) n p n+1 -p n 



At 



At 



1 ' 9(ph V,P n XP n+1 -P n ) + ^V,P n Mp n+1 -v n ) 



At I 3p 



3p 



+ ((p n+1 - p n ) 2 ) + O ((p n+1 - p n ) 2 ) 
9(ph). n n ,p n+1 -p n 



9p 



P ,P 



At 



+ 0(At) 



9 (P h ) r p n )pn) [u n . Vp n + p n v . ^ + 



ap 

9(ph) 



9p 



;p n ,p n ) [u n+1 ■ Vp n + p n V ■ u n+1 ] + O(At). 



(4.22) 



Similarly, using (4.16) we have 



u n+l . V ( p n+1 h n ) = u n+l . y^n) + Q(At) 
,n+1 



u 



r p n ,p n )Vp n + — — (p n ,p n )Vp n 



9p 



9p 



+ O(At) 



^^(p n ,p n )u- +1 -Vp n + 0(At). 
9p 



(4.23) 



Adding (]4.22[) and ( |4.23[ ) in view of ( |4.20[ ) leads to 

P 



h 



9 (ph.) 



9p 



V-u n+1 = O(At), 



or 



(|)"v.u-=0,At). 



(4.24) 
(4.25) 



With (g) ^ 0, we deduce that V • u n+1 = O(At) which ends the proof. 

The proof of the asymptotic preserving property for the fully discrete scheme follows 
the same methodology and is left to the reader. 



5 Numerical results 

In this part we provide numerical results for the second-order asymptotic preserving 
scheme, the first-order scheme being too diffusive. We first test the accuracy and the 
convergence order of the scheme on the colliding acoustic pulses test-case, then study the 
behavior of the scheme in the compressible regime with shock tubes test-cases, using the 
Euler equations. At last, we test the behavior of the scheme at low Mach number with the 



16 



well-known test-cases of the backward facing step and the lid driven cavity, modeled by 
the full Navier-Stokes equations, and the test-case of the heat-driven cavity, which uses 
the energy equation. The results are compared to the results of the Low Mach Roe scheme 
described in [T6] [T7] , using the OVAP code to run the simulations [32] . The Low Mach 
Roe scheme is an incompressible solver and has been the object of previous validation. 



5.1 Colliding acoustic pulses 



This test-case proposed in [30] consists in two acoustic pulses, a right-running pulse and 
a left-running pulse. The pulses first collide and superpose, with a maximum of pressure 
at t = 0.815s, then separate to return to their initial configuration at T = 1.63s. The 
boundary conditions are periodic. The Mach number is £ = i and we use a perfect gas 



of constant y = 1 .4, the equation of state being p 
[— L, L] is considered, with L = - 



Y V 
y-1 h 

and is discretized into 220 cells 



A one-dimensional domain 
The time-step is 



At = 1 x 1 3 s. The parameter oc in the numerical scheme is taken as a = 0. The initial 



data for this case are: 



1 x 

p(x,0) = p + ^PiU - cos(27t-)), po = 0.995, pi = 2.0 

1 x 

p(x,0) =Po+ 2 £ PiH -cos(27t-)), p = 1.0, pi=2y, 

1 x 

u(x,0) = -sign(x)u (1 - cos(27t-)), u = 2y/y. 



(5.1) 
(5.2) 
(5.3) 




(a) t=0.185s 



(b) t=1.63s 



Figure 2: Pressure profile for the colliding pulses test-case. The initial profile is in dashed 
line, and the solid line gives the result of the second-order scheme for two different times. 



The pressure profile computed by the second-order scheme is compared to the initial 
condition on fig. [2j At t = 0.815s, the pressure reaches a maximum value as the two 
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pulses are superposed. At t = 1.63s, the pulses are separated from each other again. As 
explained in [30J, weakly nonlinear acoustic effects distort the final profile as shocks are 
beginning to form in the vicinities of the locations x = ±18.5. 



Convergence tests We check that the scheme has indeed a second order convergence 
in space. We calculate the error between the solution p obtained for the pressure with 
N = 100, 200 and 400 cells with a reference solution p re f calculated with N Te f = 3200 cells. 
In order to check the spatial convergence only, the time-step is taken as At = 0.05 x Ax 2 . 

The error ||E|| L i is the discrete L 1 norm of the difference between the solution p and 
the reference solution p ref : 



II 1 = 



2_r=r \v 



ref 



Xj)l 



(5.4) 



where p(x,) is calculated by linear interpolation when x, is not a discretization point for 
the discrete solution, as p has been computed with less cells than p Tef . 

The results have been computed with oc — 1 to avoid spurious oscillations due to the 
shocks forming in the vicinities of the locations x = ±18.5. The L 1 norm of the relative 
error between the reference solution and the results for 100, 200 and 400 cells is given in 
table [Tj The logarithm of the error as a function of the logarithm of the space step Ax is 
plotted on fig. [3] We indeed have a second order convergence in space. 



Cells 


Ax 


At = 0.05 x Ax 2 


IIEIIli 


100 
200 
400 


0.44 
0.22 
0.11 


9.68 x 10^ 3 
2.42 x 10~ 3 
6.05 x 10- 4 


2.61 x 10~ 3 
7.94 x 10~ 4 
2.55 x 10~ 4 



Table 1: Colliding pulses test-case. L 1 norm of the relative error between the reference 
solution computed with 3200 cells and the numerical results for 100, 200 and 400 cells. 



5.2 Shock tube problems 

Sod shock tube This shock tube test-case has been proposed in [51]. The initial state 
is divided in a left part < x < 0.5 and a right part 0.5 < x < 1 , the initial values being 
given in table [2] 

We use a perfect gas of equation of state p = -zyj^, with a constant y = 1.4, Neumann 
boundary conditions, and a mesh with 100 cells, with Ax = 0.01 and At = 0.001 . The 
Mach number is £ = 1 and the parameter cx is zero. 

The results are represented at t = 0.2s on the left column of fig. [4] for the density, 
pressure and velocity computed by the second order scheme. The exact solution is also 
displayed. The second-order scheme shows a small overall deviation from the reference 
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_9 I i I i I i I i I i I 

-2,5 -2 -1,5 -1 -0,5 

log( dx ) 

Figure 3: Convergence test for the first order in time and second order in space scheme 
on the colliding pulses test-case. The solid line is the log of the L 1 error as a function of 
the log of Ax. The dashed line is a line of slope 2. 



solution and satisfactory shock velocities are obtained, as expected from a conservative 
scheme. 

We give the error of the solution of the second-order scheme, compared with the exact 



solution, in table. |3j The calculus of the error is given in equation (5.4). In order to 
check the spatial convergence only, the time-step is taken as At = Ax 2 . We can see that 
the presence of discontinuities reduces the space convergence order from 2 to 1 as in 



Lax shock tube This one dimensional shock tube proposed in [35] presents stronger 
shocks than in the Sod shock tube problem. The initial state is divided in a left part (I 
subscript) for —1 < x < and a right part (r subscript) for < x < 1, the initial values 
being given in table [4] 

We use a perfect gas of equation of state p = with a constant y = 1.4, Neumann 

boundary conditions, and a mesh with 100 cells, with Ax = 0.01 and At = 0.001. The 
Mach number is £ = 1 and the parameter a is zero. 





P 


h 


u 


Left 


1 


3.5 





Right 


0.1 


2.8 






Table 2: Sod shock tube. Initial conditions for the pressure, the enthalpy and the velocity. 
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Cells 


Ax 


At = Ax 2 


HElIu 


100 


0.01 


1 x 10- 4 


1.3 x 1CT 2 


200 


0.005 


2.5 x 10- 5 


6.8 x 1CT 3 


400 


0.0025 


6.25 x 1(T 6 


3.4 x 1CT 3 



Table 3: Sod shock tube. L norm of the relative error between the exact solution and 
the numerical results for 100, 200 and 400 cells. 





P 


h 


u 


Left 


3.528 


27.748 


0.698 


Right 


0.571 


3.3997 






Table 4: Lax shock tube. Initial conditions for the pressure, the enthalpy and the velocity. 



The results are shown at t = 0.25s on the right column of fig. [4] for the density, 
pressure and velocity computed by the second order scheme. The exact solution is also 
displayed. As in the Sod shock tube, the accuracy is satisfactory and shock velocities are 
accurately restored. 

We give the error of the solution of the second-order scheme, compared with the exact 
solution, in table. |5} The calculus of the error is given in equation (5.4). We can see that 
the presence of discontinuities reduces the space convergence order from 2 to 1 as in 



Cells 


Ax 


At = Ax 2 


IIEIIli 


100 
200 
400 


0.02 
0.04 
0.005 


4 x 1CT 4 
1 x TO" 4 
2.5 x 10- 5 


1.2 x 10- 2 
6.1 x 10- 3 
3.4 x 10- 3 



Table 5: Lax shock tube. L norm of the relative error between the exact solution and 
the numerical results for 100, 200 and 400 cells. 



These test-cases demonstrate the satisfactory behavior of the scheme in the compress- 
ible regime. 

5.3 Backward-facing step test-case 

The backward-facing step test-case is a two-dimensional test-case which checks the accu- 
racy of the scheme in the low Mach regime. The geometry of the step creates a region 
of low velocity where a recirculation of the fluid takes place. The size of the circulation 
region depends on the Reynolds number of the flow. This test-case has already been 
treated experimentally and numerically, for example in [2"| I5B"]. 
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(a) Density (b) Density 




(c) Pressure (d) Pressure 




x 



(e) Velocity (f) Velocity 

Figure 4: Density, pressure and velocity profiles for the Sod shock tube (left column) and 
Lax shock tube (right column). The exact solution is displayed in dashed line and the 
result of the second-order scheme (100 cells) is in solid line. 
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This case is modeled by the full Navier-Stokes equations ( 2.1[ )-(2.3). The contribution 
of the diffusive part (right-hand side terms in the equations) is added in explicit source 



terms as mentioned in part 3.1 



Figure 5: Backward-facing step - Geometrical features 



The geometry of the step is such that Li = 4m, l_2 = 1 8m, h, = 2m, where the notations 
refer to fig. [5] The Computational domain is discretized with a uniform Cartesian grid 
of step Ax = 0.2m, and the time-step is At = 5 x 10~ 4 s. The results are displayed at 
T = 20s. 

The global Mach number is e = 0.01. We take cx = 0. A perfect gas of equation of 
state p = -^y ^ and constant y = 1 .4 is used. The initial conditions are p = 1 x 1 5 Pa, 
H = 3.5 x 10 4 J/kg, u= (1,0) m/s. 

The coefficients in the diffusive terms of the Navier-Stokes equations are y = 1 .56 x 
10~ 2 m 2 /s, A = 2.7 x 10~ 2 W/m/K, C p = J/K/kg, with R = 8.315 J/mol/K and 

M. = 0.02897 kg/mol. The corresponding Reynolds number is Re ~ 75. The external 
forces are neglected (f ext = 0). 

A wall slip boundary condition (u ■ n = 0) is applied on the step and on the top and 
bottom walls. At the inlet, the velocity and enthalpy are imposed, while a Neumann 
condition is applied on the pressure. The value of the inlet velocity is u = (1,0) m/s and 
the imposed enthalpy is h. = 3.5 x 10 4 J/kg. At the outlet, only the pressure is imposed 
with a value of p ou tiet = 1 x 1 5 Pa. 

The modulus of the velocity and the streamlines computed by the second-order Asymp- 
totic Preserving scheme at t = 20s are displayed on fig. [6] for the first 10m of the channel, 
whose total length is 22m. The results are compared to the results obtained with a clas- 
sical Roe scheme and with the Low Mach Roe scheme mentioned in the introduction of 
this section. 

The second-order Asymptotic Preserving scheme gives a very satisfactory result as the 
recirculation is computed and matches the dimensions of the recirculation computed by 
the Low Mach Roe scheme. On the other hand, we can see that the Roe scheme, without 
a low Mach number treatment, is not able to capture the recirculation of the fluid. This 
test-case thus confirms that the Asymptotic Preserving scheme has a satisfactory behavior 
in the incompressible regime. 

5.4 Lid-driven cavity test-case 

The two-dimensional lid-driven cavity test-case is also a well-known problem to assess 
the adequacy of a code to the low Mach number regime (see for example [20]). The 
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(a) Second-order Asymptotic Preserving scheme 




case concerns a cubic cavity full of fluid where all the walls are immobile but one : this 
moving wall drags the neighboring fluid, which initiates a global circulation of the fluid. 
We expect a central primary recirculation and a smaller lower right eddy. 

This case is also modeled by the full Navier-Stokes equations (2.1)-(2.3), and the 
contribution of the diffusive part (right-hand side terms in the equations) is added in 



explicit source terms as mentioned in part 3.1 



The global Mach number is £ = 0.01. We take <x = 0. A perfect gas of equation of 
state p = -^y ^ and constant y = 1 .4 is used. The initial conditions are p = 1 x 1 5 Pa, 
h = 3.5 x 1 4 J/kg, u = (0, 0) m/s. 

The coefficients in the diffusive terms of the Navier-Stokes equations are v = 2.5 x 1 0~ 2 
m 2 /s, A = 2.7 x 10- 2 W/m/K, C p = J/K/kg, with R = 8.315 J/mol/K and 

M = 0.02897 kg/mol. The external forces are neglected (f ext = 0). 

The cavity is formed by the domain cu = [0,1] x [0,1], discretized with a uniform 
cartesian grid of step Ax = 1 /50. A wall slip boundary condition (u • n = 0) is applied on 
all the walls except the top wall. The top wall is moving at a speed varying continuously 
from u = at t = 0s to u = Im/s for t > Is. The computation has been run with a 
time-step of At = 2.5 x 10~ 4 , until a final time of t = 20s. 

The modulus of the velocity and the streamlines computed by the second-order Asymp- 
totic Preserving scheme at t = 20s are displayed fig. [7j The results are compared to the 
results obtained with a classical Roe scheme and with the Low Mach Roe scheme 



23 




(a) Second-order Asymptotic (b) Low Mach Roe scheme (c) Roe scheme 

Preserving scheme 




(d) Second-order (e) Low Mach Roe (f) Roe scheme 

Asymptotic Preserving scheme 

scheme 



Figure 7: Lid-driven cavity test-case - Modulus of the velocity and streamlines. 



24 



While the classical Roe scheme displays no recirculation whatsoever, the second-order 
Asymptotic Preserving scheme shows a good behavior as the circulation region is com- 
puted and is similar to the circulation region computed by the Low Mach Roe scheme. 
The primary vortex is clearly visible for the second-order Asymptotic Preserving scheme 
and the Low Mach Roe scheme on the figures displaying the streamlines, while it is not 
correctly computed by the classical Roe scheme. We can also see the lower right eddy 
expected along with the primary vortex. 



5.5 Heat-driven cavity 

The heat-driven cavity test-case consists in a two-dimensional steady-state single-phase 
laminar flow resulting from a natural convection created by the difference of temperatures 
between the two vertical walls of a cubic cavity and the gravity field. The horizontal walls 
are adiabatic walls. 

This test-case is well suited to evaluate the behavior of a numerical scheme in the 
low Mach number regime and in the presence of thermal conductivity terms and gravity 
terms. It has been studied for example in [TTJ ESI US]- This case is very interesting in our 
situation because it requires the energy equation, contrary to the two previous cases that 
could have been run with the isentropic Navier-Stokes equations. At last, the global Mach 
number resulting from the scaling of the equations is £ = 10~ 4 , which is much smaller 
than in the two previous cases. 



This case is modeled by the full Navier-Stokes equations (2.1 )-(2.3 ).A perfect gas of 
equation of state p = £ and constant y = 1 .4 is used. The dimension of the cubic 
cavity, L = 1 .528 x 10~ 3 m, is chosen so that the flow is a low Mach number flow (e = 10~ 4 ) 
; viscosity and conductivity are chosen so that the flow is laminar (low Reynolds number 
: Re ~ 37) and results from natural convection. 

The coefficients in the diffusive terms of the Navier-Stokes equations are: v = 1 .61 9 X 
1CT 6 m 2 /s, A = 2.29 x KT 3 W/m/K, C p = J/K/kg, with R = 8.315 J/mol/K and 

M. = 0.02897 kg/mol, the external force is gravity: f cxt = (0, — 9.81 )m/s 2 . 

The initial conditions are p = 1 x 10 5 Pa, h = 2.9167 x 10 5 J/kg, u = (0,0) m/s. 
A wall slip boundary condition (u • n = 0) is applied on all walls. The velocity of the 
walls is zero. The top and bottom horizontal walls are adiabatic walls : the thermal flux 
is imposed to be zero. The temperature is imposed on the right and left vertical walls : 
Ti = 283.1 5K on the left wall and T r = 263.1 5K on the right wall. 

The domain is discretized with a uniform cartesian grid of step Ax = L/40. The 
computation has been run with a CFL of 0.002, until a t = 2s and with a convergence 
criteria for the Newton method of 10~ 7 . Then the computation has been continued with 
a CFL of 0.001 and a convergence criteria of 10~ 10 until a final time of t = 3s. This case 
has been computed with the second-order Asymptotic Preserving scheme, with cc = 0. 
In order to compare the results, we also present the results of computation of the Roe 
scheme and of the Low Mach Roe scheme. We expect to find specific patterns in the 
visualization of the isocontours of the local Mach number and the temperature. 

The isocontours of the local Mach number (which is different from the Mach number e 
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(a) Second-order Asymptotic (b) Low Mach Roe scheme (c) Roe scheme 

Preserving scheme 



Figure 8: Heat-driven cavity test-case - Isocontours of the Mach number 



resulting from the scaling of the equations) is given on fig. |8} We can see that the solution 
computed by the second-order Asymptotic Preserving scheme matches the solution of the 
low Mach Roe scheme. On the other hand, the Roe scheme is not able to provide the 
correct solution and the pattern is very different from the pattern obtained with the 
scheme adapted to low Mach numbers. Let us also notice that the local Mach number 
ranges from 10~ 5 to 10~ 9 , which is very small and confirms that the case lies in the 
incompressible regime. 

6 Conclusion 

The aim of this paper was to provide an all-speed scheme for the numerical simulation of 
mixed compressible and incompressible fluid flows. The second-order discretization of the 
proposed Asymptotic Preserving scheme shows a very good behavior in both flow regimes. 
In compressible situations, we obtain good shocks properties as the scheme is conservative. 
In the low Mach number regime, the Asymptotic Preserving property provides a consistent 
discretization of the incompressible model, the divergence-free condition on the velocity 
is respected and the pressure is solved via an elliptic equation. The centered spatial 
discretization of the implicit pressure term allows the time-step to be based on the fluid 
velocity and not on the acoustic velocity. The time-step can be much larger than with an 
explicit upwind method and does not depend on the Mach number. The proposed scheme 
therefore shows a very good behavior on the weakly compressible numerical test-cases 
such as the backward-facing step and the lid-driven cavity as it provides the expected 
recirculations of the fluid, and also provides the correct solution on the heat-driven cavity 
which uses the energy equation. 

Low Mach number regimes are often encountered in multiphase mixtures. The Navier- 
Stokes equations have been used in this paper as they are very similar to the simplest 
two-phase flow model, the homogeneous equilibrium model. In future works, we intend to 
extend the scheme to more elaborate two-phase flow models as the four-equation mixture 
model and the six-equation two-fluid model. 
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First tests have been realized so far with the four-equation mixture model and a test- 
case of a water flow in a heated channel has been computed. It has confirmed the ability 
of the scheme to compute a two-phase mixture, to use a general equation of state (Water 
and Steam EOS), and to work with heat transfer terms and phase change phenomena. 
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